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Abstract 

Primordial Black Holes (PBHs) that form from the collapse of density perturbations are more 
clustered than the underlying density field. In a previous paper, we showed the constraints that 
this has on the prospects of PBH dark matter. In this paper we examine another consequence 
of this clustering: the formation of bound systems of PBHs in the early universe. These would 
hypothetically be the earliest gravitationally collapsed structures, forming when the universe is 
still radiation dominated. Depending upon the size and occupation of the clusters, PBH merging 
occurs before they would have otherwise evaporated due to Hawking evaporation. 

PACS numbers: 04.70.Bw, 97.60.Lf, 98.80.Cq 
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I. INTRODUCTION 



Primordial black holes (PBHs) are a unique probe of cosmology, general relativity, and 
quantum gravity PBHs form from the gravitational collapse of density perturbations that 
are of order unity on the scale of the cosmological horizon 

HQ. 

Measurements of the 

cosmic microwave background (CMB) anisotropy [3!] imply that density perturbations at 
the time of decoupling are much smaller (5h ~ 10~ 5 ). As such, PBH formation will be 
cosmologically negligible during and beyond this era. Less constrained are the conditions in 
the early universe before decoupling, and we cannot preclude the existence of much larger 
density contrasts which could have formed PBHs. 

PBHs would be the first gravitationally collapsed objects in the universe. As clustering is 
ubiquitous in other, observed gravitationally collapsed systems (galaxies, clusters of galaxies, 
superclusters, etc), it will be no different for PBHs. In a previous paper 4j], we derived the 
basic properties of PBH clustering and what consequences it had on the viability of PBHs as 
dark matter (DM). As PBHs are created with an isocurvature component, this constrained 
the mass and abundance of PBHs if they are to serve as DM. 

The aim of this work is to continue the analysis of PBH clustering, particularly to inves- 
tigate the creation and behavior of bound systems, or clusters, of PBHs. In Sections [III and 
III II we describe general properties of PBHs and their clustering, respectively, that will be 
used throughout the paper. The main results of the paper are presented in Section [TV] and 
conclude in Section [V] We use units throughout such that c = 1. 

II. PBH BASICS 

A black hole of mass M has a Schwarzschild radius R$ = 2GM = j^-. Throughout 
we assume that any PBHs have negligible angular momentum and electric charge, unless 
otherwise noted. It is assumed that PBHs are formed at a fraction / of the horizon mass 

M PBH = fM H , (1) 

M H (t) = M P (J-^ = (2 x 1O 5 M ) (J^j (2) 

in the radiation dominated era. An important parameter in the clustering and abundance 
of PBHs is 

v = 6 c /a rad (r H ) (3) 



which measures the height of the density peak at the threshold for PBH formation in units of 
the variance of the radiation density perturbation a ra d smoothed over the horizon radius r# 
when that perturbation enters the horizon. We will examine the consequences of clustering 
for a range of values of u, f and M#, as there is considerable variation in their respective 
values. The horizon fraction / depends upon details of the radiation perturbation profile 
and equation of state (see a discussion in \4\). The horizon mass Mh and v are determined 
primarily by the power spectrum peak location and value, respectively, the details of which 
will depend upon the exact PBH formation model invoked. Possibilities include: models 
where PBHs are produced at a particular scale in pre/reheating during/after multi-field 
(hybrid) inflation (s 9 ] ; models with a peak at a specific mass/length scale due to a drop 
in pressure during a phase transtion (see [? ? ] for examples relating to the QCD phase 
transition); generic inflationary models with blue spectra [? ? ]. In all of these models, the 
PBH mass function is strongly peaked, and so for simplicity we take a monochromatic mass 
function. 

The fluctuation of a density field (either for radiation or PBHs) is defined in terms of the 
mean density, 

S-'-?. (4, 

PBH production is exponentially suppressed (e 2"); exact expressions for the initial den- 
sity depend upon v and assumptions as to how the overdense region is defined. Assuming 
the peaks are Gaussian, the probability of a given point in the smoothed (over the horizon 
scale) density field falling above the PBH formation threshold is 

- 1 - ■ > -1/2 / <P 

'Be V ' V 2cr Ld( r H), 

which will roughly be equal to the PBH number density at formation across many horizon- 
sized regions. In the limit that v ^> 1, we obtain an approximate expression by taking the 
upper limit to infinity, 

P = erfc ( -^=) a (J^vA e-» 2 l\ (6) 



ft = I (27ra 2 rad (r H )) exp - 2 - - dS, (5) 
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More precise calcuations (assuming, e.g., that the given point is not just above threshold but 



also is a maximum) for a power law (of slope n) radiation spectrum 10 
prefactors for the above exponential. We define the function N*(y) by 



11] give different 



d = N^v)e-^l\ (7) 



with values for different models being 



J-v 1 erfc approximation , 

NJv) = { , \\ /2 \. (8) 



The initial PBH density is then 

N4v)e-» 2 / 2 

npBH = t^ = r? (9) 

Vh Vh 

in terms of the Hubble volume at formation 

V H = (2^ 2 R 3 H . (10) 

Due to quantum effects Jjj], a BH of mass M will emit particles as a blackbody with 
temperature T h given by 

T h (M) = — L_ = w 10 - _1 eV (11) 

" l j SttGM 8ttM U ^lgj 6V - 1 j 

As the temperature is inversely proportional to the mass, this is unobservable for a one solar 
mass BH (Th(M & ) w 62 nK), but cannot be neglected in the mass range of PBHs. This 
emission also corresponds to a mass loss for the PBH, 

M = -L h = -o* 8B Tt{4*B*) = (12) 

where a* SB is the effective Stefan-Boltzmann constant and is proportional to the effective 
number of relativistic degrees of freedom in the emitted particles. PBHs therefore have a 
finite liftime, after which they would have emitted their entire rest mass, given by 

(10" s) — (13) 



3«(M ) v J \lgJ 

As the lifetime scales with M 3 , there is a threshold mass above which holes will not have 
evaporated by the present day (to)- This threshold mass M* is given by 



M* fa (4 x 10 14 g 



n(.U.) \f /„ 



(14) 



QM x 10 25 g 3 /s/ V4.4 x 10 7 s, 

Given the uncertainties in a and to, a threshold mass of M* ~ 10 15 g is typically quoted in 
the literature. 

A large enough abundance of PBHs with M ~ M* will produce a number of observable 
effects through their evaporation in the current day. They would contribute to cosmic rays 



galactic center 



13 1, the 7-ray background Q, 511 keV emission due to positron annihilation in the 



16| or be the cause of short duration gamma ray bursts 17|, |l8] . Observations 



(or the lack thereof) of PBHs evaporating today depend critically upon not only the number 
density of PBHs present today npsH{to), but also upon how clustered they are within within 
the galaxy. Assuming an isothermal halo model, the effective number density is C n PBH{to) 



13 



15| and ranges from 10 5 — 10 7 . 



where £ is the local density enhancement factor 

PBHs with M < M* would have evaporated by the present day. The main mechanism 
for "observing" PBHs in cosmology is through their Hawking radiation. In the absence of 
a direct detection, the main utility of PBHs is to set limits of PBH abundance at various 
limes given a non-detection. Though, PBHs have also been invoked to explain baryogenesis 



191 ] . reionization 20|] and as a solution to the magnetic monopole problem 21]. 

Evaporating PBHs have their most dramatic effect during the era of BBN, where Hawking 



radiation can alter the entropy per baryon and light element abundances [22|? ? ]. Therefore, 
the success of BBN implies an upper limit to the number of PBHs evaporating at that time. 
Combining Equations (12]), (JTJ) and ( ITBl gives the relation 

f 3 M P / ^ 3 



r(t) 



3a \tp 



(15) 



the lifetime r of a PBH created at a time t. What this allows one to do is use information 
from a "late epoch" (time r) to examine conditions at an "early epoch" (time (<t). In 
the above example, r ~ tsBNi an d the limits on initial PBH abundance from BBN imply 
(3 < 1(T 16 for M PBH between 10 9 g and 10 15 g. 

This relation depends critically upon the PBH only losing mass through evaporation, 
and not gaining mass in any way (accretion or merging). Should this not be the case, the 
lifetime r is no longer given by the initial PBH mass, and the link between late epoch and 
early epoch is broken. Instead, the energy in PBHs that would have evaporated away can 
now linger for longer periods of time. Since r oc M 3 , the merging of two equal mass BHs 
will result in a BH with a lifetime 8 times as long. If this merging can continue, then there 
is a greater chance of PBHs produced in the early universe still existing today. 
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III. PBH CLUSTERING 
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Here we briefly summarize the results of |4|. Defining the radiation and PBH two point 
correlation functions £(r) and power spectra P(k) under the assumption of spherical sym- 
metry: 

£rad(r) = (S rad (x)S rad (x + r)) (16) 

£,pim{r) = (5pbh(x)S PB h{x + r)) (17) 

Using a Gaussian window function Wk(krn), we will work with smoothed versions of the 
correlation functions over the horizon size r#: 

Ud(r) = ^J dkk 2 P rad {k)*-^^\W k (kr H )\ 2 (18) 

Zpbh(t) = / dkk 2 P PBH (kf-^^\W k (kr H )\ 2 . (19) 
It is useful to redefine the radiation field correlation function in a normalized fashion 

[ ) ~ ~ oUth) (20J 

so that w(0) = 1. 

As PBHs form at the peaks of the radiation density field (where a radiation perturbation 
5 > 5c), the PBH correlation function is itself a function of the radiation correlation function. 
In general, for Gaussian perturbations, the peaks of a density field are more clustered than 



the density field itself; this was first illustrated by Kaiser [23( comparing rich clusters of 
galaxies (Abell clusters) to just galaxies. This will be the case for PBHs too; PBHs will be 
more clustered than the underlying radiation field {^pbh » £rad)- 

As originally shown in 24|, [25|], the two-point correlation function in the limit of large v 



becomes [4j 

1 + 6W0 =exp(u 2 w(r)). (21) 

Note that the radiation density field must possess large (order unity) fluctations at the 
horizon scale at which PBHs are formed; this implies that the PBH "fluctuations" must 
be even larger (exponentially so). Their evolution is thus intrinsically nonlinear. As the 
universe expands and nearby PBHs enter the same Hubble volume, they immediately detach 
from the general expansion and collapse to form bound clusters. Their dynamics will be 
determined by the size of the cluster and typical separation within the cluster. As PBHs 



form in (initially) separate horizon volumes, the horizon distance si the smallest possible 
distance scale over which PBHs could be correlated. As we will show later, the closest pairs 
of PBHs will form withing a few horizon distances from each other. 

The exact details of the PBH clustering depend upon the underlying density field, rep- 
resented here by the exact form of w(r). Given $,PBH{ r ) we can transform to find the 
PBH power spectrum Ppbh{^)- Because PBHs are discrete objects, the limit of the power 
spectrum as k — > is a constant and related to the number density of PBHs. The power 
spectrum for a group of N uniformly randomly distributed objects over a volume V 3> Vh is 
1/N = {ripBuV)' 1 = Our Pp_B_ff(0) 7^ indicating that the PBHs are distributed as 
clusters of objects with mean occupation number N c = Ppbh{0)P- Due to the exponential 
enhancement of the correlation function, N c scales exponentially with v as well. The exact 
relation will depend upon the assumed form of w(r), and thus the assumed power spectrum 
Prad(k)- In {4], we assumed a Gaussian spike and obtained N c ~ N*(v)e v ' !//4 . 

IV. PBH CLUSTER EVOLUTION 
A. Literature 



PBH clusters were first explored in Freese et al. [26|], where PBHs in the mass range 
M € [10 15 g, 10 33 g] form clusters after matter-radiation equality. The collapse of baryons 
onto these clusters creates explosions that act as the seeds of large scale structure in the 
Ostriker-Cowie [27j model. 

As shown in the previous section, PBH perturbations enter the horizon with a very large 
amplitude (<Jpbh ~ e^ 2//2 ). It is therefore no longer valid to treat their evolution using linear 
perturbation theory, as one is able to do for other forms of CDM. Instead, we examine the 
sub-horizon evolution of the PBH population as an N-body problem. As noted earlier, using 
numerical simulations to examine PBH formation is difficult given the small numbers of them 
that form. This is more so true should one want to examine their subsequent non-linear 
evolution after creation. 

We instead appeal to previous work done in the context of DM halo formation and N- 
body simulations to examine the PBH population behavior. Being non-relativistic, PBHs 
will cluster hierarchically (just as CDM); creating bound systems that get incorporated into 
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larger ones. The internal dynamics of these systems are determined solely by gravitational 
clustering, analagous to other gravitationally bound systems such as star clusters and galax- 
ies. For this, we are aided by the work done in the context of studying more massive black 



holes in globular clusters 



28 



and galaxies 



291 ] . In those cases, gravitational interactions 



tend to either produce bound pairs or ejections, rather than BH coalescence [30|. What 
occurs in the case of PBHs depends upon how many (N) form in a "PBH cluster" and what 
their initial separations D are. 

We shall discuss the clustering in a hierarchical manner as well, first discussing the small 
scale dynamics (PBH binary formation and evolution), then the larger scale (cluster relax- 
ation, ejection, and evaporation). 



B. PBH Binaries 



PBHs will form at rest with respect to the background expansion. A PBH perturbation 
entering the horizon implies that PBHs are able move appreciably in response to the gravity 
of other PBHs. The cluster as a whole begins with no kinetic energy and, properly defined, 
negative potential energy, for negative total energy. Starting at rest also implies no net 
angular momentum (either orbital or in BH spin) for the cluster, and further no net angular 
momentum for any subset of the cluster. The first physical process to consider is direct 
PBH merging due to a head-on collision between two PBHs in a cluster (to create a PBH 
with sum of the two initial PBHs' masses). For the simplest case of N = 2, this type 
of merging is automatic, as there are no significant tidal forces to keep the PBHs from 
colliding 1 . For any larger cluster (N > 2), the situation becomes more complicated. In this 
case, any pair (or subset) of PBHs will have non- vanishing orbital angular momentum with 
respect to any other member of the cluster, and tidal forces suppress the head-on merging 
of PBHs. Rather, close encounters between PBHs are more prone instead to create PBH 
binary systems through a dissipative process via emission of gravitational waves. We will 
address more carefully the conditions for binary formation in Sec JIV Gl 

Thus, we expect binary PBH formation to be common within the cluster, and we examine 
the evolution of PBH binaries. That PBHs would form binaries in the RD era was for 

1 Insignificant tidal forces would arise from density perturbations in the background radiation field. While 
a "nearby" perturbation might be just below threshold for PBH creation, it will decay upon horizon entry. 
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studied in the context of MACHO PBHs of around a solar mass [3 lM34| due to the expected 
gravitational wave emission. There are a number of other physical processes other than 
gravity that may affect the PBH binaries. To estimate the magnitude of each effect, we first 
assume that two PBHs of identical mass M are in a circular orbit of radius D (a separation 
of 2D). In the initial limit that the PBHs are well separated compared to their Schwarzschild 
radii (Rs = 2GM <C D) and also moving non-relativistically (v << 1), the Newtonian force 
is 

_GM^_Ml(RW 
9rav ~ AD 2 ~ 1Q \ D ) [ ] 

The orbital timescale ( "period" ) is therefore 



/ D \ ^^ 2 / R \ ~~ 

t orbit = AttD ( —) = 4V2ttD -4- 



(gm) = ^ d [t) (23) 

We next examine three mechanisms that afect the orbital motion of a PBH binary: 
background radiation drag, Hawking radiation pressure, and gravitational wave emission. 

C. Background Radiation Drag 

PBHs moving through the background will accrete radiation, and this will have a com- 
pensating effective drag force on the PBHs. Considering the radiation to be a perfect fluid 
with sound speed c s = cj y/3, the accretion rate for an individual PBH will be: 



M bondi = p r irR acc ^v 2 + c 2 s , (24) 

where v is the PBH physical (not comoving) velocity with respect to the rest frame of the 
radiation and R acc is the accretion radius, defined as 



Race — -5— 2 ( 25 ) 

Due to the radiation sound-speed being so high, the accretion radius can be approximated 
by its non-relativistic limit of R acc = 3Rs- 

Thus, the drag force F drag ~ M hondi v ~ Ma, or 
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giving the time-scale for drag 

M 



t, 



drag 



V 

a 



M 



Bondi 



1 _Pr_ 

3 3 / 4 • 47r I Mp 



M" 1 = — 1^1 — 



x/3 /M 2 \ 1 



1J 



p 

Pr 



(27) 



Note that this is also the mass-doubling time for the PBH due to accretion. In the radiation 
dominated era we can rewrite this as 



^drag 



\/Z /8tt 
"18 V~3 



1 

Rs 



16 / t 



t. 



(28) 



Since Rs <tn = 2t at creation, td rag > t and therefore radiation drag is negligible 



35|. 



D. Hawking Radiation Pressure 

PBHs will emit Hawking radiation with a temperature T^. This leads to a finite PBH 
lifetime, as previously discussed. Further, the Hawking Radiation pressure: 

Ph = a SB T^ (29) 

where a* SB is the (dimensionless) modified Stefan-Boltzmann constant for Hawking radiation. 
This leads to a radiation force of 

And a timescale for orbital disruption of 



•—Hnss-J^lt 1 1 (31) 



note this is different from the PBH lifetime r; 

/ 2 5 7T 3 



T 



M 2 P R 3 S (32) 



So that 

U i- 2 21/2 / D \ 3/2 

The conclusion being that the PBH will completely evaporate well before the evaporated 
radiation pressure can disrupt the binary. 
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E. Gravitational Radiation 



Finally, For close binaries of PBHs, the loss of energy due to gravitational radiation may 
become important. We estimate the size of this effect as follows. Consider two identical 
mass PBHs separated from their center of mass by a distance D in a circular orbit. The 
emitted power due to gravitational waves is 

P = 19^1. (34) 
5 R 5 K J 

Given a gravitational binding energy U = GM 2 /(2D), the time-scale for inspiralling due to 
gravitational wave losses is 

U _ 5 / D \ 3 
tspiral ~P~4 [gMJ { ' 

This can roughly be modelled as an additional drag force; 

f - - - JS- (**V /2 m) 

wave ~ v ~ 5 • 25/2 \ d J 

Gravitational wave emission will cause the binary orbit to decay and accelerate the PBH 
merging process. Comparing to the orbital timescale, 

twave _ 5^/2 / D \ (37) 

t orbit 47T \RsJ 

For D > Rs, t wave t orbit and, as expected, many orbits occur before inspiraling. 



F. Conditions for Binary Merger 

Therefore, the only mechanism that will affect the PBH binaries will be gravitational 
wave emission. The concern now is whether the PBHs can evaporate before the inspiral is 
complete. Comparing the relevant timescales, 

, '"- lh '» GB '~ (0.166) r^V^V (38) 



r ~ 2% 3 R e s v ' \ M J \R S J 

If t spiral /t < 1, then the PBH binary will merge before evaporation takes place. This limit 
can be turned into a maximum value of the initial separation D for which merging will 
occur. To gauge what importance this effect will have, we need to know more about the 
distribution of initial separations. 
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First, define two quantities. From the definition of S(x) (equation H]) , p(x) = p(l + S(x)). 
Consider the quantity 

(p(x)p(x + f))=p 2 (l + a^)), (39) 
which can be shown by direct substitution and noting that (S(x)) = (5(x+r)) = 0. Assuming 
isotropy and a monochromatic mass function, Ppbh = np BB M, and 

(n PBH {x)n PBH {x + f)) = n 2 {l + £(r)). (40) 

The mean initial separation we define as D = fipgjj. This, being the first moment of the 
density distribution, tells us nothing about clustering. Define instead D = (np BH (0))~ 1 ^, 
which has the dimensions of length. Using equation [21] in the limit that v » 1 and r 
(as we are interested in the clustering at short distances), we find 

(n 2 p BH (0)) = n 2 PBH e» 2 (41) 

Using Equation (JS)) and the definition of np BB , 

D = V2^[N,(iy)]-^ 3 e u2/(i R H (42) 

so that 

D = {n 2 PBH {^))-^ = De-» 2 " = sFH[N*{v)]- 1 '*R h (43) 

While the number density is exponentially suppressed with u, there is a cancellation 
computing this clustering length so that it is roughly constant with v. As expected from the 
peak-background split model, the scale that sets the clustering is the horizon size at PBH 
formation. 

Taking D = D in Equation 1381 for / = 1 and N*(v) = J~v~ x , we obtain 



^ = 2.66zW|^) 2 , (44) 



t \M H , 

meaning, provided M H > M P and with only slight v dependence, inspiral typically occurs 
well before Hawking evaporation completes. In Figure [1] we plot the conditions for merging 
before evaporation for the range of hole masses expected to evaporate by the current day. 
We see that, assuming binary formation occurs, the PBH clustering implies certain merging 
over most of this range, even relaxing the assumption on /. 

Let us talk a bit more about this result. For PBHs with M < 10 15 g, this implies that if 
the formation of PBH binaries is efficient, the lifetime of the PBH population as a whole is 
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15 



nu=1 



D_bar 




10 3 10 u 10° 10 lu 10 

M (g) 



FIG. 1: Conditions for binary merger before evaporation. Df, ar is the mean separation assuming 
a uniform (non-clustered) PBH distribution. D tw id = (f v 1 / ? 'f~ 1 Rs is the "clustering length" 
scale. Solid curves are for / = 1, dashed curves are for / = 10~ 35 . 

increased, as r oc M 3 . Given hierchical clustering, where a new generation of PBH binaries 
are formed from the first generation, this process can theoretically proceed until all smaller 
PBHs are bound up into larger ones. 

This would have a profound impact upon cosmology - as PBHs that should have evapo- 
rated by now would still be present in the universe. There are two caveats in this scenario. 
First, if PBH binary formation is inefficient, then most PBHs will still evaporate "when they 
should" and the merged ones have little impact. Second, it assumes that the PBH binaries 
are not disrupted by other means; such as other PBHs. We address this point next. 

G. PBH cluster occupation 

Having shown that gravitational wave emission can cause PBHs to merge before they 
evaporate, the creation and destruction of PBH binaries in a PBH cluster will determine 
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whether merging is a common enough even to impact the PBH mass function as a whole. 

We define a PBH cluster as a group of PBHs that are within causal contact (sub-horizon) 
and are gravitationally bound. We showed in [4j that, under the assumption that v > 4, 
PBHs are distributed at large distances as clusters of objects with mean occupation number 

N c ~ N*(v)e v2/A . (45) 

The question then remains: what range of parameters (M , v) are physically interesting, 
given the absence of evidence for PBHS in our current universe? In [4], we constrained 
these according to the criteria that a) PBHs do not induce an early matter-dominated phase 
of the universe, and b) PBHs, as dark matter, do not produce a too large isocurvature 
perturbation. These give a parameter space where the allowed v decreases as Mo increases 
(see their Figures 1-3). 

We can then look at two different categories of PBHs, distinguished by mass: those that 
would and would not have evaporated by the current day, with the boundary at M ~ 10 15 g. 
The subsequent evolution of this first category (the lightest PBHs) is interesting in that the 
evaporation lifetime r oc M 3 - if merging occurs, this will increase the PBH lifetime, perhaps 
above the "evaporation" boundary. For Mh < 10 15 g, this gives a range of v ^ 7 — 10 for 
/ = 10 -3 ' 5 and v ^ 8 — 11 for / = 1. For Mh > 10 15 g there is no worry of Hawking 
evaporation, but these have a smaller lower bound on v and thus would form (on average) 
smaller clusters. 

The evolution of bound PBH clusters will depend upon their internal density profile. The 
average cluster profile is just (ppbh(j)) = ppbh (1 + (^pbh^)))- Since a randomly chosen 
point is likely to be located at a peak, it can be shown that[? ]: 

< W0> = r^TnTW(r) « &bh{t). (46) 

The cluster is highly peaked, due to the enhancement of the PBH correlation function. 
The exact details of the shape of the cluster will depend upon choice of initial radiation 
power spectrum P ra d{k), window function Wk(kr) and higher order statistics (we are only 
considering the 2-point function). Nonetheless, while PBHs can form no closer than a 
horizon distance apart, their typical separation within the cluster D cannot be too much 
larger than this distance due to the exponential enhancement of the correlation function for 
small r. Due to the high central density (below we compute the concentration parameter c), 
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this results in the PBHs being approximately "close-packed" inside the core. Outside the 
core, the density falls off faster than the underlying radiation perturbation, as £pBH(r) = 
exp (z/ 2 u>(r)) — 1. Far from the core, the cluster profile falls off no slower than the underlying 
radiation perturbation profile: as r — > oo, £p_B#(r) ~ v 2 w(r). The steepest profile comes 
from assuming P ra d oc ^Dirac^ — ^*) ( a delta function a given scale which we can take 
to be the horizon scale at formation); there w(r) = sm ^ fc * r ^ and asymptotically w{r) oc r _1 
as r — > oo. In this extreme limit, normalizing the cluster profile (ppbh^)) to the mean 
objects in a cluster N c puts virtually all PBHs within k' 1 - this would be unphysical give 
this is roughly the event horizon size, but motivates our later assumption of putting all of 
the PBHs within the core. A more realistic power spectrum flattens this profile; taking 
P Ta d oc k n gives w{r) oc r~( n+3 ) asymptotically [? ]. 

The quasi-equilibrium state of an N-body system of gravitating point masses can be 
broken down into two major components: a high density central core and a low-density 
encompassing halo 



30]. Due to the enhanced clustering, the PBH clusters considered here 
already begin with a similar, centrally concentrated profile. It is convenient to discuss 
separately activity within the core (small radius, or r — > 0) and in the halo (large radius, or 
r oo) of the cluster. The core concentration c is then 

c = ESSH = 1 + !^i^ PBH{0) - = o pbh - e^ 2 . (47) 

PKalo l + gff$6W00) 1 

This alone is a clue that our cluster is destined for collapse; an isothermal sphere with 
c > 708.61 is unstable to collapse (known as the gravothermal catastrophe) [? ]. This 
happens for v > 3.62, which is in the parameter space we were considering anyway. Further, 
as the core begins to contract and PBHs begin to approach relativistic speeds, the relativistic 
instability [? ] accelerates this collapse. To confirm this, we examine in more detail the 
processes involved in this collapse process. 

As our clusters begin "pre-collapsed" in a sense, there is expected to be a great difference 
between core radius and halo radius when it comes to estimating the extent of the clusters. 
Due to its high density, it is within the core that most interactions take place [? ? ], and 
so it is within the core we are more interested anyway. We can make some estimates of size 
using our expressions for N c and D, assuming the majority of PBHs form within the core. 
As N c oc V C ore oc -R?ore' the cluster radius R cor e oc N^ 3 . In the core, where the PBHs are 
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most concentrated, each PBH is separated by a typical distance of D, then 

1/3 



Rcore ~ N^D = (iV,(z/)e^ 2 / 4 ) x s/2^(N*(u))- 1/3 R H = V^e^l l2 R H . (48) 

Thus, the core size is independent of our choice of N*(i/), but increases with increasing 
v, as expected. Taking a range of 2 < v < 10, this gives a range of 3.5Rh < Rcore < IO^Rh- 
As the higher values of v are only allowed for PBHs that form earlier, these large cores will 
consist of smaller mass PBHs (and a smaller horizon size), while smaller cores will have 
larger mass PBHs. 

We can similarly make an estimate for the initial size of the surrounding halo. Recall that 

— —1/3 

the mean separation of PBHs is D = n P ^ H . As we are now grouping PBHs into clusters of 
size N c , the density of clusters of PBHs is smaller an the density of PBHs themselves, thus 
the mean separation between clusters D c must be greater than the mean separation between 
PBHs themselves: 

D c = (n PBH /N c y 1/3 = N X J 3 D ~ (iV,(z/)e^/ 4 ) 1/3 x (N*(v)y 1/3 e» 2 / 6 R H = V^e^^Rn 

(49) 

This sets an upper limit on the halo radius, imagining the PBH clusters close-packing in 
space. As the evolution of the cluster is driven primarily by core activity, we will focus on 
that from here on out. 

Take a PBH cluster of iV initial PBHs of mass M with scale radius R c = sRs, where 
Rs = 2GM is the initial Schwarzschild radius. The dynamical timescale (roughly the cluster 
crossing time) is 

^ffrr^ s3nRs - (50) 

As a reference, the binary period and inspiral timescales are 

t orhit = 4TrV2d 3/2 R s , (51) 
t spiral = WD =10d 4 R s , (52) 

where d = D/R s . 

The first process we need to consider for large N-body systems is that of relaxation. This 
is the process by which a cluster achieves equilibrium (virialization) through the combined 



effect of two-body scatterings. The relaxation timescale is given by 

/ 0.14iV 



36? 



rel 



\n(0AN) 
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t c . (53) 



Gravitational wave emission during 2-body scatterings will accelerate cluster relaxation (so 
that the expression for t re i is an upper limit), in addition allowing for binary mergers. Under 
the assumption of violent relaxation, however, the cluster virializes in only a few dynamical 
times. After virialization, the PBHs will have a velocity dispersion 

. 2 , GM C N 

<" > ~ = s (54) 

The simplest mechanism for PBH mergings is direct collisions; the timescale being 

t coll ~ 0.8 ln(0.4iV) (j^q) Uei (55) 

where = (4u 2 ) -1 . This process is initially negligible compared to the relaxation time for 
small (N < 100) clusters. 

Successive scatterings can give a single object enough energy so that it can escape from 
the cluster entirely (v > v escape ), while the cluster shrinks in size. This mass loss gradually 
leaves to the "evaporation" of the cluster; where most of the bodies are ejected to infinity, 
leaving behind only a hard binary system (in the absence of gravitational wave emission) 
or a central black hole (where gravitational wave emission has induced orbital decay and 
merging of the hard binary). The timescale for this evaporation is 

tevap ~ 300t re ;. (56) 

The cluster shrinking is accompanied by core collapse where the innermost portions of 
the cluster accrete more and more of the mass of the cluster, which could result in runaway 
growth unless halted by some mechanisms. Numerical simulations [? ? ] show that this 
process begins within 10-20 relaxation times - much smaller than the cluster evaporation 
timescale. 



H. Binary Formation 

We now turn to binary formation in a cluster. This is a critical point, because a popula- 
tion of binaries in the core could eject other PBHs from the core while themselves contract- 
ing in their orbit (binary hardening). This, effecively, cools the core, and could possibly 
halt core collapse (see [36| for a discussion). This is the case when PBHs are moving 
non-relativistically (in the Newtonian regime). Once relativistic effects are considered, in 
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particular gravitational wave emission, we see that this is not the case, and binaries are not 
enough to arrest core collapse. 

Since the PBHs form at rest with respect to the cosmic expansion, there are no binaries 
present at the "birth" of the cluster 2 . In the Newtonian regime (cluster members mov- 
ing non-relativistically), the only avenue for binary formation is through 3-body exchange 
interactions only, with a formation timescale of 



As cluster "evaporation" occurs in a few hundred relaxation times, binary formation through 
3-body interactions is negligible for large N systems (N ^ 100). 

Binaries are only formed through 2-body interactions when there is some energy dissipa- 
tion involved during (what would have been) a scattering event. For stars, this is through 
dissipation in the stellar atmospheres. For black holes, this happens due to gravitational 
wave emission during the scattering, which heretofore we have not considered. This has 
been studied in detail for compact clusters of black holes by a number of authors: in [? ? 
], analytic esimates and Fokker-Planck simulations are presented, while N-body simulations 
are presented in [? ? ]. While their simulations were for stellar-mass black holes, the results 
are mass independent. There, 2-body formation is dominant over 3-body formation, with 
the ratio of timescales (from [? ]): 



At the outset of cluster contraction, v < 1 and 2-body binary formation is dominant for 
iV > 300 clusters. Only as the core evolves can v — > 1 (i.e., becomes relativistic) and N 
shrink enough for 3-body formation to become important. Despite the enhanced formation 
in the core, the inspiral time t spira i < t re \ is smaller than the cluster relaxation time, so that 
these binaries are not effective at heating the core and halting collapse (noted in [? ]). 

Just what mass fraction f core of the initial cluster ends up in a central black hole is not 
well known. Numerical calculations by [? ? ] find f core ~ 0.06 — 0.1, though the simulation 
begins to break down at that point, meaning this is a lower limit. Of those PBHs not 
captured in the core (either ejected from the cluster or remaining in the halo), they will 

2 Known as primordial binaries, where primordial is being used in a similar context. 



t 3 ~ 10iV 2 In Nt re i 



(57) 



*3 




(58) 
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either Hawking evaporate (if they are light enough) or remain until the universe becomes 
matter-dominated and they are bound up within galaxies. The possibility remains that some 
residual number of PBHs from this process might survive until the current day, providing 
an observational test of this scenario. 



V. CONCLUSIONS 



A consequence of the PBH clustering, developed in the previous section, is the merging 
of PBHs into more massive, longer lived PBHs during the radiation dominated era. This 



implies that PBHs are strong candidates to be t 
distinct from other models in the literature 
which are of two types. First 



37 



37 



he "seed" BHs that form SMBHs. This is 
where PBHs serve as SMBH "seeds", 



39j | , PBHs formed from a blue spectrum of perturbations 



undergo hierarchical merging along with DM halos in the MD regime. It was found, however, 
that it is still difficult to account for the observed BH mass without invoking some additional 
accretion source; in this case the accretion of a cosmological quintessence field. Quintessence 
fields are types of scalar fields that arise in theories of dark energy, and unlike a cosmological 



constant, BHs are able to accrete energy from a time-dependent scalar field 
type of PBH seed theory [38 



40|. The second 



assumes that the ~ 1OOOM seed BHs are PBHs formed at just 
that right mass scale. This requires a deviation from scale-invariance in the power spectrum 
very close to the era of BBN, which is highly constrained. 



The advantage of our model is that PBH formation occurs much earlier than in 
so that the power spectrum isn't as constrained. Further, PBH merging takes place in 
the RD epoch, so there is more time for PBHs to merge, and possibly no need to invoke 



A more detailed study of this model is 



an additional accretion mechanism as in [37J, 
planned to determine whether this is the case. 

PBH merging in clusters dramatically changes the limits on initial PBH abundance f3. 
The rate of merging is sensitive to the initial conditions of the cluster. We have considered 
an idealized scenario, with PBH formation happening at a single mass scale and at a single 
time. Physically, one would need to account for PBH creation across a span of times, and 
include in the cluster dynamics the effects of a spectrum of masses. N-body simulations 
would need to be carried out to examine this further. 

In addition to providing the seeds of SMBHs, this PBH merging scenario we have dis- 



19 



cussed has other predictions. One prediction is more gravitational wave emission than 
originally assumed for a uniform PBH population. This is due to the increased probability 
of PBH binary formation and emission from N > 2 bound states, specifically in the core of 
a PBH cluster. This would alter recent predictions of gravitational wave spectra from direct 
graviton emission from PBHs 
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